Cacoa analysis walkthrough

Author: Katarina Dragicevic (PhD student), Biotech Research and Innovation Centre, University of Copenhagen

This notebook is about how the Cacoa pipeline was ran on different Conos objects (P15/P60 excitatory/inhibitory) The same pipeline was used for all objects and we show here on 1 example how it was done.

1. Load Conos object and annotation file. Construct the sample groups objects. Load color palettes.

#we changed some of the annotation names after we ran Cacoa already on the objects. This is why we use this recvec to recode the names for the final plots to make them match those in the paper

recvec <- (setNames(c(‘Sst_Crhr2’, ‘Sst_Crh’, ‘Sst_Nr2f2_Cdh9’, ‘Sst_Hpse’, ‘Sst_Nr2f2_Glra3’, ‘Pvalb_Kcnk10’, ‘Pvalb_Gpc6’, ‘Pvalb_Sema5a’, ‘Pvalb_Gpr149’, ‘Pvalb_Tacr1’, ‘Pvalb_Calb1’, ‘Sst_Necab1’, ‘Vip_Chat’, ‘Pvalb_Htr2c’, ‘Vip_Col15a1’, ‘Id2_Car10’, ‘Id2_Car2’, ‘Id2_Ano3’, ‘Sst_Calb2’, ‘Sst_Piezo2’, ‘Sst_Nr2f2_Necab1’, ‘Sst_Chodl’, ‘Pvalb_Vipr2’, ‘Lamp5_Lhx6’, ‘Vip_Thsd7b’, ‘Lamp5_Lsp1’, ‘Vip_Sgcz’, ‘Vip_Gpc3’, ‘Lamp5_Dock5’, ‘Pvalb_Gabrg1’, ‘Lamp5_Ndst4’), ci60$cell.groups %>% unique %>% as.character()))

recvec

               Sst_Crhr2                  Sst_Crh           Sst_Nr2f2_Cdh9                 Sst_Hpse 
             "Sst_Crhr2"                "Sst_Crh"         "Sst_Nr2f2_Cdh9"               "Sst_Hpse" 
         Sst_Nr2f2_Glra3      Pvalb_Sema3e_Kcnk10               Pvalb_Gpc6        Pvalb_Reln_Sema5a 
       "Sst_Nr2f2_Glra3"           "Pvalb_Kcnk10"             "Pvalb_Gpc6"           "Pvalb_Sema5a" 
            Pvalb_Gpr149     Pvalb_Tacr1_Il1rapl2       Pvalb_Calb1_Necab1               Sst_Necab1 
          "Pvalb_Gpr149"            "Pvalb_Tacr1"            "Pvalb_Calb1"             "Sst_Necab1" 
                Vip_Chat Pvalb_Calb1_Necab1_Htr2c              Vip_Col15a1                Id2_Car10 
              "Vip_Chat"            "Pvalb_Htr2c"            "Vip_Col15a1"              "Id2_Car10" 
                Id2_Car2                 Id2_Ano3                Sst_Calb2         Sst_Mme_Fam114a1 
              "Id2_Car2"               "Id2_Ano3"              "Sst_Calb2"             "Sst_Piezo2" 
        Sst_Nr2f2_Necab1                Sst_Chodl              Pvalb_Vipr2               Lamp5_Lhx6 
      "Sst_Nr2f2_Necab1"              "Sst_Chodl"            "Pvalb_Vipr2"             "Lamp5_Lhx6" 
       Vip_Thsd7b_Prss12               Lamp5_Lsp1                 Vip_Sgcz                 Vip_Gpc3 
            "Vip_Thsd7b"             "Lamp5_Lsp1"               "Vip_Sgcz"               "Vip_Gpc3" 
             Lamp5_Dock5             Pvalb_Gabrg1         Lamp5_Ntn1_Ndst4 
           "Lamp5_Dock5"           "Pvalb_Gabrg1"            "Lamp5_Ndst4" 

#library(cacoa)

#library(qs)

#con <- qread(“inhibitory_p60.qs”, nthreads = 5) #conos object

#anno <- qread(“anno_inhibitory_p60.qs”, nthreads = 1) #annotation

#sg.pal <- qread(“sg.pal_inhibitory_p60.qs”, nthreads = 1) #palette for sample groups

#cg.pal <- qread(“cg.pal_inhibitory_p60.qs”, nthreads = 1) #palette for cell groups

#sg <- qread(“sg_inhibitory_p60.qs”, nthreads = 1) #sample groups

#how does sg look like?
sg

           1014_WT           1099_15q P60_AAM1033_15q_b3 P60_AAM1035_15q_b3 P60_AAM1037_15q_b3  P60_AAM11a_15q_b1 
                wt                15q                15q                15q                15q                15q 
  P60_AAM11a_wt_b1  P60_AAM11b_15q_b1   P60_AAM11b_wt_b1     P60_AAM3_wt_b1     P60_AAM8_wt_b1  P60_AAM978_15q_b3 
                wt                15q                 wt                 wt                 wt                15q 
  P60_AAMT2_15q_b2   P60_AAMT3_15q_b2    P60_AAMW1_wt_b2    P60_AAMW2_wt_b2 
               15q                15q                 wt                 wt 
Levels: 15q wt

#how does sg.pal look like?
sg.pal

          15q            wt 
"dodgerblue2"       "black" 

#how does cg.pal look like?
cg.pal

        Id2_Ano3        Id2_Car10         Id2_Car2      Lamp5_Dock5       Lamp5_Lhx6       Lamp5_Lsp1      Lamp5_Ndst4 
       "#e07ba9"        "#a9e07b"        "#857be0"        "#a5e07b"        "#7be096"        "#e0c97b"        "#e07bc4" 
     Pvalb_Calb1     Pvalb_Gabrg1       Pvalb_Gpc6     Pvalb_Gpr149      Pvalb_Htr2c     Pvalb_Kcnk10     Pvalb_Sema5a 
       "#7b92e0"        "#a0e07b"        "#7be0d4"        "#e07b87"        "#afe07b"        "#ac7be0"        "#e07b94" 
     Pvalb_Tacr1      Pvalb_Vipr2        Sst_Calb2        Sst_Chodl          Sst_Crh        Sst_Crhr2         Sst_Hpse 
       "#e0b67b"        "#e0997b"        "#8fe07b"        "#cc7be0"        "#7bc1e0"        "#b5e07b"        "#e07b82" 
      Sst_Necab1   Sst_Nr2f2_Cdh9  Sst_Nr2f2_Glra3 Sst_Nr2f2_Necab1       Sst_Piezo2         Vip_Chat      Vip_Col15a1 
       "#a67be0"        "#e0c27b"        "#cbe07b"        "#7be0c2"        "#7bb3e0"        "#e07bc3"        "#7b91e0" 
        Vip_Gpc3         Vip_Sgcz       Vip_Thsd7b 
       "#7be08e"        "#7b9ae0"        "#e0c17b" 

UMAP

#cao$plotEmbedding(color.by = “cell.groups”, font.size = 3)

UMAP colored by sample

#cao$plotEmbedding(color.by = “sample”, font.size = 0)

Cell group proportions

#cao$plotCellGroupSizes()

Cell loadings

#cao$estimateCellLoadings()

p2 %>% ggplot( aes( x = values, y = ind, color = ind)) +
  geom_violin(scale = "width", 
                #fill = "#ededed", 
                color = NA, aes(fill = ind), alpha = 0.2) +
  theme_light() + theme(legend.position = "none") + scale_fill_manual(values = palcomp) +
  geom_vline(xintercept = 0, linewidth = 0.1) +
  new_scale_fill() + 
  stat_summary(fun.data = iqr, geom = "errorbar", show.legend = F, 
               aes(color = ind), 
               stroke = 0.5, alpha = 0.9,color = "grey30", width = 0.25) + 
  stat_summary(fun.data = iqr, geom = "point", show.legend = F,
               aes(fill = ind), 
               shape = 23, stroke = 0.5, color = "grey30", alpha = 0.9, size = 2.5) +  
  scale_fill_manual(values = palcomp)

Warning: Ignoring unknown parameters: `stroke`

Cell Density

#cao$estimateCellDensity(method = “graph”)

#cao$estimateDiffCellDensity(type = “wilcox”, n.cores = 200)

cao$plotDiffCellDensity(type = "wilcox", min.z = 0, adjust.pvalues = T,  contours = c("Sst_Nr2f2_Glra3", "Pvalb_Vipr2"), contour.conf = "20%")

  |                                                                                                    |   0%, ETA NA
  |================================================                                                 |  50%, ETA 00:02
  |=============================================================================================| 100%, Elapsed 00:03

Expression shifts

#cao$estimateExpressionShiftMagnitudes(top.n.genes = 300)

p3 %>% ggplot( aes( y = value, x = Type, color = Type)) +
  geom_violin(scale = "width", 
                #fill = "#ededed", 
                color = NA, aes(fill = Type), alpha = 0.2, 
              trim = F) +
  theme_light() + theme(legend.position = "none", axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + scale_color_manual(values = pal2)+ scale_fill_manual(values = pal2) +
  geom_hline(yintercept = 0, linewidth = 0.1) +
  stat_summary(fun.data = iqr, geom = "errorbar", show.legend = F, aes(color = Type), shape = 23, stroke = 0.5, alpha = 0.9,color = "grey30", width = 0.25) + 
  stat_summary(fun.data = iqr, geom = "point", show.legend = F,aes(fill = Type), shape = 23, stroke = 0.5, color = "grey30", alpha = 0.9) 

Warning: Ignoring unknown parameters: `shape` and `stroke`

Cluster free expression shifts

#cao$estimateClusterFreeExpressionShifts(n.cores = 200, gene.selection = “z”, normalize.both = T, n.top.genes = 300)

DEG analysis

#cao$estimateDEPerCellType(n.cores = 100, name = “edgeR”, n.cells.subsample = 5, resampling.method = “bootstrap”, n.resamplings = 20)

Gene ontology

#cao$estimateOntology(type = “GSEA”, name = “GSEA”, org.db =org.Mm.eg.db, verbose = T, n.cores = 200, de.name = “edgeR2”)

#cao$estimateOntology(type = "GSEA", name = "GSEA", org.db =org.Mm.eg.db, verbose = T, n.cores = 200, de.name = "edgeR2")
cao$plotOntologyHeatmap(name = "GSEA", genes = "all")

Loading required package: DOSE

DOSE v3.24.2  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/

If you use DOSE in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609